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COMPUTER SUBROUTINES FOR THE ESTIMATION OF NUCLEAR REACTION 
EFFECTS IN PROTON-TISSUE-DOSE CALCULATIONS 

John W. Wilson and G. S. Khandelwal* 

Langley Research Center 

SUMMARY 

Calculational methods for estimation of dose from external proton exposure of arbi- 
trary convex bodies are briefly reviewed and all the necessary information for the estima- 
tion of dose in soft tissue is presented. Special emphasis is placed on retaining the effects 
of nuclear reaction, especially in relation to the dose equivalent. Computer subroutines 
to evaluate all of the relevant functions are discussed. Nuclear reaction contributions for 
standard space radiations are in most cases found to be significant. Many of the existing 
computer programs for estimating dose in which nuclear reaction effects are neglected . 
can be readily converted to include nuclear reaction effects by use of the subroutines 
described herein. 


INTRODUCTION 

When an object is exposed to external radiation, the dose field within the object is 
a complicated function of the character of the external radiation, the shape of the object 
(including orientation), and the object's material composition. Calculation of dose within 
an object involves solution of the appropriate Boltzmann transport equation in which the 
external radiation source imposes boundary conditions on the solution. Although general 
purpose computer programs exist for making such estimates (ref. 1), they are seldom 
used in practice when the object is bounded by a complicated surface, as is the human 
body, for example. Instead, calculations are usually made for simple geometric shapes 
from which inferences are then made for more general geometries, and the resultant 
errors are uncertain. 

In the case of external proton radiation, such as that encountered near high-energy 
accelerators, in space, and in high-altitude aircraft, it is the inclusion of nuclear reaction 
effects which presents the major hurdles in an accurate calculation. It was found in ref- 
erence 2 that dose estimation could be greatly simplified and still include the effects of 
nuclear reaction. Furthermore, it was shown that the method of reference 2, when in 
error, was always conservative. Required for such calculations is a knowledge of the 
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transition of protons in semi -infinite slab geometry, which is the simplest geometry for 
existing transport computer programs. Indeed, almost everything that is known about the 
dose in humans from external proton radiation is inferred from calculations with slab 
geometry (ref. 3). 

In the present report, a general method for estimation of dose in arbitrary convex 
geometry in terms of dose conversion factors in slab geometry is briefly discussed. 

These dose conversion factors for protons in tissue are then represented using buildup 
factors. A parametric form for the buildup factors is presented. The values for the 
parameters are derived from Monte Carlo calculations of various authors. All the infor- 
mation necessary to estimate dose and dose equivalent for proton irradiation of convex 
objects of arbitrary shape is contained herein. The advantage of the method is that exist- 
ing proton dose calculations in which nuclear reactions are neglected can be directly con- 
verted by substitution of the dose conversion factors presented herein. 

SYMBOLS 

aj buildup -factor parameter, i = 1, 2, 3, 4 

D(x) dose at point x, rad or rem 

E proton energy, MeV 

E 0 proton energy parameter, MeV 


E r reduced proton energy, MeV 

F(z,E) proton dose buildup factor, dimensionless 

P(E) nuclear survival probability in tissue 

p(E) particle rigidity, MV/c 

p Q particle rigidity parameter, MV/c 


Qp(S) 


quality factor, dimensionless 


r(E) proton range in tissue 
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e(z) energy of proton with range z in tissue, MeV 

ct(E) proton macroscopic cross section in tissue, cm - * 

t(E) proton total optical thickness, dimensionless 

<p(£2,E) proton differential fluence spectrum, protons/cm^-MeV-sr 

<p 0 proton differential fluence spectrum parameter, protons/cm^-MeV-sr 

fi unit vector in direction of proton motion, dimensionless 

Abbreviations: 

GCR galactic cosmic radiation 

SCR solar cosmic radiation 

A prime indicates a variable of integration. 
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THEORY 


In passing through tissue, energetic protons interact mostly through ionization of 
atomic constituents by the transfer of small amounts of momentum to orbital electrons. 
Although the nuclear reactions are far less numerous, their effects are magnified because 
of the large momentum transferred to the nuclear particles and to the struck nucleus it- 
self. Unlike the secondary electrons formed through atomic ionization by interaction with 
the primary protons, the resulting radiations of nuclear reaction are nearly all heavily 
ionizing and generally have large biological effectiveness. Many of the secondary parti- 
cles of nuclear reactions are sufficiently energetic to promote similar nuclear reactions 
and thus cause a buildup of secondary radiations. The description of such processes 
requires solution of the transport equation. The approximate solution for the transition 
of protons in 30-cm-thick slabs of soft tissue for fixed incident energies are presented 
in references 4 to 11. The results of such calculations are dose conversion factors for 
relating the primary monoenergetic proton fluence to dose or dose equivalent as a func- 
tion of position in a tissue slab. 

Whenever the radiation is spatially uniform, the dose at any point x in a convex 
object may be calculated according to reference 2 by 


D(x ) = C f R n |z x (ft),E| <*>(«,E) dftdE (1) 

J 0 J ft 

where R n (z,E) is the dose at depth z for normal incident protons of energy E on a 
tissue slab, 0(ft,E) is a differential proton fluence along direction ft, and z x (ft) is the 
distance from the boundary along ft to the point x . It has been shown that equation (1) 
always overestimates the dose but gives an accurate estimate when the ratio of the proton 
beam divergence due to nuclear reaction to the body radius of curvature is small. Equa- 
tion (1) is a practical prescription for introducing nuclear reaction effects into calculations 
of dose in geometrically complex objects such as the human body. The main requirement 
is that the dose conversion factors for a tissue slab be adequately known for a broad 
range of energies and depths. 

Available information on conversion factors is for discrete energies from 100 MeV 
to 1 TeV in rather broad energy steps and for depths from 0 to 30 cm in semi-infinite 
slabs of tissue (refs. 4, 5, 8, and 9). The nuclear reaction data used for high-energy 
nucleons are usually based on Monte Carlo estimates (refs. 12 to 14) with low-energy neu- 
tron reaction data taken from experimental observation. The quality factor defined by the 
International Commission on Radiobiological Protection (ref. 15) is used for protons. The 
quality factor for heavier fragments and the recoiling nuclei is arbitrarily set to 20, which 
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is considered conservative, although the average quality factor obtained by calculation is 
comparable to estimates obtained through observations made in nuclear emulsion (ref. 16). 

To fully utilize equation (1), the fluence-to-dose conversion factors for normal inci- 
dent protons on a tissue slab must be known for all energies and depths. A parametriza- 
tion of the conversion factors was introduced by Wilson and Khandelwal (ref. 2) which 
allowed reliable interpolation and extrapolation from known values. In the following sec- 
tion, a refinement and extension of that work will be discussed. 

Fluence-to-Dose Conversion Factors 

The conversion factor R n (z,E) is composed of two terms representing dose due to 
the primary beam protons and the dose due to secondary particles produced in nuclear re- 
action. Thus, 

R n (z,E) = R p (z,E) + R s (z,E) (2) 

where the primary dose equivalent conversion factor is given by 

R p (z,E) = P(E) Q F [s(E r )] (3) 

in equation (3) is calculated by using 


(4a) 

where 

z average tissue atomic number 

A average tissue atomic weight 

I adjusted ionization potential 

m e electron mass 

e electron charge, C 


The linear energy transfer (LET) denoted by S(E) 
Bethe's formula above 243.8 keV as given by 


S(E) - 


47rN^e^z 

r 

f 2m e v2 

A 

V 2 

m e v 2 A > 

%\ 

1 

C 2 


V 

L \ c2/ ! 

J 
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V 


proton velocity, cm/sec 


c velocity of light, cm/sec 

Na Avogadro's number 


At proton energies below 243.8 keV, the LET is calculated by the empirical expression 

S(E) = E°- 303 (2517 - 6283E) (4b) 

which approximately accounts for electron capture and the inner shell corrections in soft 
tissue. The proton range in soft tissue is given by 

r(E) = \ (5) 

J 0 S(E’) 

with the reduced energy in equation (3) given by 

E r = e[r(E) - z] (6) 


where e(z) is the inverse function of r(E). The total nuclear survival probability for a 
proton of energy E is given by 


P(E) = exp 


_r E 

J o. 


cr(E’) 


dE' 

S(E') 


( 7 ) 


where the macroscopic cross section cr(E) for tissue as calculated by Bertini is given 
in reference 17. The proton total optical thickness given by 

t(E) = f cr(E') -^1 (8) 

•>0 S(E’) 

is tabulated in table 1 for purposes of numerical interpolation. In the case of conversion 
factors for absorbed dose, Rp(z,E) is taken as 

Rpfe.E) = P(E) ^1 (9) 


6 


Buildup Factors 


The representation of the conversion factors is simplified (see ref. 2) by rewriting 
equation (2) as 


R n (z, E ) - 


l + 


Rg(z,E) 


Rp(z,E) 


Rp(z>E) 


or 


R n (z,E) = F(z,E) Rp(z,E) 


( 10 ) 


where F(z,E) is recognized as the proton dose buildup factor. The main advantage for 
introducing the buildup factor into equation (10) is that unlike R n (z,E), the buildup factor 
is a smoothly varying function of energy at all depths in the slab and can be approximated 
by the simple function 

F(z,E) = (aj + a 2 z + a 3 z 2 ) exp (-a 4 z) (11) 

where the parameters a^ are energy dependent. The aj coefficients are found by fit- 
ting equation (11) to the values of the buildup factors as estimated from the Monte Carlo 
calculations of proton conversion factors. The resulting coefficients are shown in table 2. 
The coefficients for 100, 200, and 300 MeV protons were obtained by using the Monte Carlo 
data of reference 4. The values at 400, 730, 1500, and 3000 MeV were obtained from the 
results of Alsmiller, Armstrong, and Coleman (ref. 8). The 10 GeV entry was obtained 
from the calculations of Armstrong and Chandler (ref. 9). The values that are footnoted 
in table 2 were obtained by interpolating between data points or smoothly extrapolating 
to unit buildup factors at proton energies near the Coulomb barrier for tissue nuclei 
(=12 MeV). The resulting buildup factors are shown in figures 1 and 2 and are compared 
with the Monte Carlo results; the error bars were determined by drawing smooth limit- 
ing curves to bracket the Monte Carlo values following the general functional dependence. 
These uncertainty limits should, therefore, be interpreted as limits of approximately two 
standard deviations, rather than the one standard deviation usually used in expressing un- 
certainty limits. 


CONVERSION FACTOR COMPUTER CODE 

To utilize equation (1) in a specific problem requires values for the conversion fac- 
tor R n (z,E) over the range of interest. Formulas for these factors were presented in 
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the previous section. A computer code has been generated to return values of R n (z,E) 
for arbitrary depth z and energy E. This code is listed in the appendix and is described 
briefly here. There are six main functions to be generated relating to LET, range -energy 
relations, quality factor, and the functions relating to nuclear reaction effects given as 
nuclear survival probability and buildup factor. 

The functions relating to ionization by the primary beam are generated by the func- 
tion subroutine RTIS. Tables for r(E) and S(E) are generated on the first call to RTIS. 
Subsequent intermediate values are found by numerical interpolation above 10 KeV. A 
simplified approximation based on equation (4b) is used at lower energies. The function 
e(z) is found by numerical inversion of r(E). 

The quality factor is approximated by 

Q f (S) ~ 0.06S 0 - 8 (12) 

for S greater than 35 MeV/cm and set to unity for smaller LET. 

t . . 

The values shown in table 1 for the optical thickness are generated in the function 
subroutine PN and stored in an array for numerical interpolation; the nuclear survival 
probability is calculated using equation (7). 

The coefficients for calculating the buildup factors are generated by subroutine 
ANTER as a function of energy by interpolating between the values shown in table 2. 

The conversion factors are generated by subroutine RESP by supplying parameters 
z and E, which represent distance in centimeters of tissue and proton energy in mega- 
electron volts; respectively. The returned values of the conversion factors have units of 
rads (or rems) per proton per centimeter squared. 

SAMPLE CALCULATIONS 

To illustrate the usage of the buildup factors described here, calculations of the dose 
in slab geometry for normal incident protons with spectra typical of the space environment 
have been made; calculations were also made neglecting nuclear reaction effects. The 
percentage contribution to the absorbed dose and dose equivalent due to nuclear reactions 
are shown in figures 3 and 4, respectively. The spectra indicated by GCR in the figures 
represent galactic cosmic radiation with the spectrum given by 

^GCR^ = (!3) 
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The spectra denoted by the parameter p Q represent solar cosmic ray spectra given as 


0SCr( e ) = ex P 


-p(E) 

Po 


(14) 


with the rigidity given as 


p(E) = I 


E^E + 2mp) 


1/2 


(15) 


where q is the proton charge and nip is the proton mass. The value p Q = 100 MV/c 
corresponds to an intermediate -energy solar event and p 0 = 400 MV/c corresponds to a 
high-energy solar event. The curve denoted by E 0 = 100 MeV represents the energetic 
inner belt protons with a spectrum 

0(E) = 0 O eX p||T^ ( 16 ) 

It is clear from the figures that dose estimates for galactic cosmic rays and high- 
energy solar cosmic rays cannot be accurately calculated without proper account for 
nuclear reactions. This is especially true for estimates of the dose equivalent. Although 
reasonable estimates (±10 percent) of low and intermediate solar cosmic ray absorbed 
doses are expected, the dose equivalent estimates must include nuclear reaction effects. 
Marginally good estimates of absorbed dose for inner belt protons can be made by neglect- 
ing nuclear reactions, but dose equivalent estimates require inclusion of nuclear reaction 
effects. 


CONCLUDING REMARKS 

A set of computer subroutines have been developed to estimate dose and dose equiv- 
alent in tissue for incident protons, and the use of these subroutines in dose calculations 
in complex geometry has been discussed. It was found that numerically the effects of 
nuclear reactions generally cannot be neglected in the calculation of doses due to space 
proton irradiation. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, Va. 23665 
May 3, 1976 
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APPENDIX 


PROGRAM LISTING FOR CONVERSION FACTOR CALCULATION 

The computer subroutines given in this appendix were developed for the present cal- 
culations, except for MGAUSS and IUNI which were taken from the mathematical subroutine 
library of the Langley Computer Complex. 
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APPENDIX 


subroutine wesp cen. x ,rad,rem > 

C THIS SUBROUTINE GENERATES VALUES FOR THE SLAB CONVERSION FACTORS 

C FOR VALUES OF PROTON ENERGY EN <M4*V) AND DEPTH IN THE SLAB X (CM) 

REAL C < b ) 
ener=en 
FX = X 

CALL ANTER (ENER *c « r ) 

RRES = R-Ex 
ENFRP=ET IS (RRES ) 

I F ( FNERP ) 3 A *33* 3a 

33 CONTINUE 
R AD = 0 • 

rem=o. 

RETURN 

34 CONTINUE 

CALL APROB < EX *KNFR . PROti ) 

CALL atoppienerp.stopp ) 

2 CALL AF<ST0PP*GALF) 

22 PES=PRO0*STOPP*QALF 
PRINT 1 .C 

C0R£0= < C ( 1 ) + X* (C (2 ) +X»C ( 3 ) ) ) *EXP ( -X*C ( 4 > > 

COREC- ( C < 5 ) + X* (C ( 6 ) + X*C ( 7 ) ) > *EXP <-X*C (B ) ) 

IF <COREQ.LT • 1 . ) COREO= 1 • 

IF (COREC.LT . 1 • ) COREC* 1. 

PRINT I ♦EV*X«EN£RP« , STOPP * OALF * COZf-.O ' CO^EC 

1 FORMAT < 2X , 8F 12*3) 

REM=PES*COREO *1.6E"8 

PES=PROB+STOPP 

RAD=PES*COREC*l .6E-8 

return 

END 


SUBROUTINE ANTER(ENER*OR> 

C THIS SUBROUTINE GENERATES THE VALUES OF THE PARAMETERS 

C OF THE ANALYTIC FITS OF THE MONTE CARLO RESULTS 

REAL C ( 8 ) , A ( 1 2 « 8 ) , £ ( 12) 

LOGICAL FALS 

DATA E/30. , 60 • , 1 00. , I 50. ,200 . , 300 • , 400 . , 730 • , 1 200 ..1500., 3000. , 
1 1 0000. / 

DATA A/ 1 .0. I .2, 1 .4, l .5. I .6, 1 • 70, 1 .90,3.40,4.32,4.60.5.35,6.20. 

2 0 .0 , O. O . • 02 , .07 , .09 , • 1 1 . . 1 3 , • 1 56 . . 1 67 , . 1 70 ♦ . I 90 , . 280 , 

30.0. .0. 0.0. 0.0 .0.0. 0.0. 0.0. .00035. .00145. .0025. .0030. .0035, 

40. C.. 013.. 030, .^285 , .040, .033 , .0228 , • 0 1 50 , . 0 1 3 , . 0 1 2 . • 0 1 0 , • C 1 0 , 

51. 0. 1. 0.1. 1.1. 12. 1.15. 1.2. 1.24. 1.4. 1.67.1.8,2.. 2. 3, 

60 .0, . 01 . .04 0, ,06, .0 62. . 069 , .071 , .09 , .094 , .095, . 1 O* • 1 1 . 

70.0. 0.0.0. 0.0.0. 0.0.0 .0.0.0, , OuC 1 • . 0C080 • .CO 1 5 , • 002 . ,00205, 

80.0. .01 . .026. .031 ..032. .026. .0226. .015. .0122. .012. .01 ..01, 

DATA FALS/.T./ 

DATA IPT/-1/ 

R=RT I S ( ENER ) 

IF(FALS) GO TO 10 
1 CONTINUE 

ELOG=ALOG(ENER) 

CALL 1UNI ( 12* 12, E* 8* A. 2, FLOG , C * I PT , 1 ERR ) 

RETURN 

10 CONTINUE 

DO 1 1 1=1.1? 

E C I ) =ALOG (E( I ) ) 

1 1 CONT INUE 
FALS=.F . 

GO TO 1 
t Nn 
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SUBROUTINE AF (STOPPtOALF ) 

C THIS SUBROUTINE COMPUTES THE QUALITY FACTOR AS A FUNCTION OF 

C LINEAR ENERGY TRANSFER 

IF (STOPP-35. ) 1 I ♦ 1 1 . 1 2 
1 \ OAI tr = i . 

RETURN 

12 GALF=.06*STOPP**.8 
RETURN 
END 


SUBROUTINE APROrf ( EX . E ♦ PROB ) 

C THIS SUBROUTINE GENERATES VALUES FOR THE NUCLEAR SURVIVAL PROBABILITY 
C OF A PROTON OF ENERGY E (MEV) AFTER TRAVELING A DISTANCE EX (CM) IN TISSUE 

RRES=RT ISCE ) -EX 
PRO0=O. 

IF (RRES.LE.O. ) RETURN 
ENEW=ET I S ( RRES ) 

PROB =PN ( E ) /PN ( ENEW ) 

RETURN 

END 


FUNCTION PN ( E > 

C PN gives probability that proton TRAVELS FULL range WITHOUT 

C being absorbed 

EXTERNAL fox 
LOGICAL TRU 
REAL R ( 30 ) . E T ( 30 ) 

DATA ET/O. . 1 u. .2B. .50. , 1 00 • • 1 50 . . 200 • *250. ,300 • ,350. .400. . 500. » 

I 70 0. , 9U0 • . 1 I UO • * 1 300 . ♦ 1 500* ♦ 1 7CQ • . 2000* *2200 • » 24 00 • . 2600 • * 2B00. « 
23000 • . 4 uO 0 • « 5 JO'-' • ♦6O00« *7000* . B500. . 10000./ 

data tru/.t./ 

DATA IPT/-1/ 

IF(TRU) GO TO 10 
1 1 1 ER=F 

CALL IUNI (30.30.ET. 1 .R.2.ER«BYRD. IPT. IERR) 

PN=EXP< -BYRD ) 

RETURN 
10 TRU=.F. 

R ( 1 ) = 0 . 

DO 1 1=2.30 

EU=ET ( I ) 

G = ET ( I - 1 ) 

CALL MGAUSS ( G.EU.04 . ANS.FOX.F. I ) 

R ( I > =R ( I - 1 ) +ANS 
1 CONTINUE 
PRINT 19 

19 FORMAT (Z//.25X , *PN GRID*//) 

PRINT 119 

119 FORMAT < 1 OX • *F VALUES FOR GRID*//) 

PRINT 226. ET 

226 FORMAT (2X.flEI6.6) 

PRINT 227 

227 FORMAT <//.lOX.*R VALUES FOR GRID*) 

PRINT 226. R 

GO TO 111 
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SUBROUTINE FOX(X.F) 


ENFR=X 

3 call AS igmcener .s igma ) 
call ATOPPIRNER. stopp ) 
F=s I GMA/STOPP 
2 RETURN 
END 


SUBROUTINE ASIGM(ENER*SIGMA ) 

THIS SUBROUTINE GENERATES VALUES OF TOTAL NCNEL AST I C MACROSCOPIC 

CROSS SECTION (CM**2/G) IN TISSUE AS A FUNCTION OF PROTON ENERGY ENE«<NEV) 

REAL EN ( A3 ) *C«OS< A3 ) 

DATA EN/25.32,29.06* 3* • 1 6 • 39 • 86 , 44 . 65 » SO • O I . 60 • 1 9 . 70 • 24 . 79 • A 7 , 89. 9 
1 1 . 100.8, 117.9.139.3* 1 56. 3. 175*3. I 86 . 6 . 202 • o . 266 . 1 . 304.7,375.2,407. 

27.471. 6 «bO 7.1 .574.5*61 1.4.67B.3.714 .5.776.4.809.3.870.4.916.3. 1 007 

3. . 1 1 29. . 1 406. . 1 736. . 2024 . . 23 1 8. . 307 1 . . 3408 • . 3943. . 5000 . . 8000 • . 

4 1 0000. / 


oata 

CROS/ 2.614 * 2 . 360 , 

2 . 153.1 

• 905 . l 

. 387. 1 

. 757 . ; 

1 . 621 . 

I .526 . 

I .451 . 

I • 

I 379 . 

l . 327. 1 .261 • 1 . 211 « 

. 1 . 187,1 

. 164,1 

[. 152.1 

1 . 141 , 

1 . 097 . 

1 . 087 , 

1 . 100 . 

1 • 

2 1 36 « 

1 . 1 99 . 1 . 21 2*1 . 266 - 

* 1 . 293* 1 

. 350 , i 

1 • 379 , ; 

1 .424 , 

1 . 440 , 

1 .471 . 

1 . 4 78 . 

1 . 

3504 . 

1 . 477. 1 . 480 • 1 . 483 . 

* 1 . 485 . 1 

. 487 , 1 

l . 475 . I 

l .461 . 

1 . 463 . 

1 . 46 * 1 

. 450 . 



4 1 .452/ 

DATA 1 PT/-1 / 

11 E=ENER 

IF <ENER.LT. 25.32 J £N£R=25.32 

CALL I UNI <43.43.EN. 1 .CROS . 2 . ENER . CPOSS . 1PT . I ERR) 

SIGMA= < CROSS/ 10O # > 

enfr=e 

RETURN 

END 


FUNCT ION R T I S < E ) 

C THIS SUBROUTINE GENERATES THE RANGE-ENERGY RELATIONS AND LET FOR 

C PROTONS IN TISSUE 

EXT CRN AL ATOE 

REAL ET <57 ) .RT (57 ) . ST (57) 

LOGICAL FALSE 
data false/.t./ 

DATA NP/57 / 

DATA ET/.01 . .02. .03. .04. .05. .06* .07. .08. .00. .1 ..2. .3. .4. .5. 

1 .6. . 7. .8. .9. 1 . .2. .3. *4. . 5 . .6.. ^ • .3. *9. . 10..80..30* ,40. .50. * 

260. . 70. .80. .90. . luO. .150. *200. *300. .4 JO. .500. . 60Q. .700- * 

3800. . 900. . 1 000 • . 1 50 o • . 2CCC • * 2500. .3000. .4000. . 50 CO • * 6COO . . 
47000. • 8500. . 1 OOuO • / 

N= 1 

lF(FALSE) GO TO lO 
12 CONTINUE 

RT [S=E#* .697/(251 7.* .697) 

IF (F.LT. .01 ) RETURN 
A = ALOG < E ) 

DO I IE=2.NP 

IF < A.LT .ET < IE ) ) GO TO 2 

1 CONTINUE 

2 I = IF 

SLOPE = <RT < I l-RTU-1 ) ) / (ET ( I )-E r < I - 1 ) ) 

WAL = RT ( l - 1 ) + SLOPE * ( A-E T ( I - 1 ) > 

RT I S-EXP ( RAL ) 

RETURN 
ENTRY ST IS 


13 



APPENDIX 


N = 2 

1F(FALSE)G0 TO 10 

13 CONTINUE 

RTIS 2 E**. 303* C 2b 1 7.-6283. *E ) 

IF (E.CT . .01 > RETURN 
AsALOGCE ) 

DO 3 I E =2 *NP 

IF < A.LT .ET< IE > > GO TO 4 

3 CONTINUE 

4 I = IE 

slope = cst < i >- st ( 1 - 1 > j/<et< i j-et < i - i > ) 
sal = ST( I-l ) + SLOPE* ( A-ET( 1-1 ) ) 
RTIS=EXP(SAl ) 

RETURN 
ENTRY ETIS 
N = 3 

IF (FALSE )G0 TO 10 

14 CONTINUE 

RTIS=<2517.*.697*E >**1 .4 34 72 
IF(E*LT..Ol) RETURN 
R=ALOG(E ) 

DO 5 IQ=2*NP 
IF (R.LT.RT( I R ) ) GO TO 6 

5 CONTINUE 

6 I = IR 

SLOPE 2 (ET< I >-ET( I-l ) )/ <RT< I )-RT ( I-i > ) 

EAl =ET ( I-l J+SLOPE* (R-RT ( I-l ) ) 
RTIS=EXP(EAL ) 

RETURN 

10 CONTINUE 
RT ( 1 ) =0 • 

ST ( 1 ) =0 . 

M = 06 

no PI I =2 ♦ NP 

CALL A T OPP C E T < I ) . ST ( I )) 

CALL MGAUSS(ET < I-l ) «FTT( I ) • M « ANS « A T OE * F • 1 ) 
21 RT ( I ) =RT (I-l ) + ANS 
RIRST=RT (2 ) 

E I RST 2 E T ( 2 ) 

DO 1 1 ' X = 2.NP 

ET ( IX )=ALOG(ET ( IX ) ) 

RT < IX )=ALOG(RT ( IX ) ) 

11 ST ( IX )=ALOG(ST ( IX > ) 

FAL5P=.F . 

GO TO ( 1 2* 1 3* 1 4 )N 
END 


SUBROUTINE A T OE ( E * F ) 

CALL ATOPP(E*S) 

F = 1 ./S 
RETURN 
END 


SUBROUTINE ATOPP<ENER.STOPP ) 

this subroutine computes the stopping power for proton in tissue 

I F (ENER.GT . .2436 ) GO TO 2 

ST OPP= (251 7. -6263.+ENER )*ENER**.303 

RETURN 

2 ZETA=ENER/938. 21 1 

BE TAS=( (ZETA*(ZETA+2. ) )/( (ZETA+l . )**2 ) ) 

WBE=1 . 02220 l£6*eETAS/( 1 .-BETAS) 

FBET=ALOG ( WBE ) -bet AS 

STOPP= .307261 48* ( -2 • 2378342+ . 529726*FBET ) /BETAS 

RETURN 

END 
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subpout ine mgauss ( a.b«n.sum,func .fofx. number > 


DIMENSION U(5) .R<5) * SUM ( i , ,fOFXI 1 ) 
DO I LL=! *NUMPFR 
1 SUM(LL)=0.0 

IFIA.EO.B) RETURN 
U ( 1 ) = . 4256628205091 Ha 
UC2 ) = .283302302935376 
U ( 3 ) = • 1 602952 1 5850468 
U<4 )= .06746631 6655508 
U < 5 ) = .0 1 304 673574 1 4 1 4 
R ( 1 ) = . 1 47762 1 1 2357376 
R ( 2 > = • 1 3463335965499 
R<3 ) =. 1 09643181257991 
R<4 1 =.074725674575290 
R (5 » = .0 333356721 54 34 4 
F IN£=N 

OELT A=F |N£/<B-A ) 


no 7 K = 1 » K' 

XI =K-l 

F 1 NE = A + X I /DELTA 
00 2 11=1.5 
UU = U ( I 1 >/DELTA+F INE 
CALL FUNC (UU.FOFX) 

DO 2 JOYBOY= 1 .NUMBER 

2 SUM( JOYBOY >=«< I 1 1 *F0FX(30YB0Y ) + SUM ( JOYBO Y ) 
00 3 JJ = 1 .5 

UU=( 1 .0-U< JJ ) J/DELTA+F INE 
CALL FUNC (UU.FOFX) 

DO 3 NN=l. NUMBER 

3 SUM (NN > =R( JJ > *F OF X ( NN ) + SUM ( NN ) 


C 

C* 

C* 

C* 

c * 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 
c* 

C* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 


DO 7 I JK = 1 ♦ NUMBER 
7 SUM ( I JK ) =SUM < I JK ) /DELTA 
RETURN 
END 

SUBROUTINE IUNI <NMAX.N.X,NTAB«Y. lOROeR'XU, YO. IPT . IERR ) 


PURPOSE9 


USE9 


DAPflMPTrRS9 


ntab 


SUBROUTINE 1 UN I USES FIRST OR SECOND ORDER 
lagRangian inteqpolat ion to estimate the values 

OF A SET OF>UNCTIONS AT A OO I NT XO • I UN I 

USES ONE INDEPENDENT VARIABLE TABLE AND A DEPENDENT 

VARIABLE TABLE FOR EACH FUNCTION TO BE EVALUATED. 

The ROUTINE ACCEPTS the INDEPENDENT VARIABLES SPACED 
AT EQUAL OR UNEQUAL INTERVALS. EACH DEPENDENT 
VARIABLE TABLE MUST CONTAIN FUNCTION VALUES CORRES- 
PONDING TO EACH X(l) IN THE INDEPENDENT VARIABLE 

table, the estimated values ARE returned IN the yo 

ARRAY WITH THE N-TH VALUE OF THE ARRAY HOLDING THE 
VALUE OF THE N-TH FUNCTION VALUE EVALUATED AT XO. 


CALL JUNI ( NMAX . N. X. NTAB « Y. i ORDER. XO. YO . IPT. I F RR ) 


THE MAXIMUM NUMBER OF POINTS IN THE INDEPENDENT” 

variable array. 

THE ACTUAL NUMBER OF POINTS IN THE INDEPENDENT 
ARRAy.WHERE N *LE« nmax. 

A ONE-D I MENS I ONAL ARRAY. DIMENSIONED (NMAx) IN THE 
CALLING PROGRAM, WHICH CONTAINS THE INDEPENDENT 
VARIABLES. THESE VALUES MUST RE STRICTLY MONOTONJC. 

THE NUMBER OF DEPENDENT VARIABLE TABLES 

A TWO-DIMENSIONAL ARRAY DIMENSIONED { NMA X • NT A0 ) IN 
THE CALLING PROGRAM. EACH COLUMN OF THE ARRAY 
CONTAINS A DEPENDENT VARIABLE TABLE 


I UN 1001 0 

* JUN 10020 

* I UN I 0030 

* I UN IC040 

* I UN I 0050 

* I UN 10060 

* I UN I 0070 

* I. UN 10080 

* 1 UN I 0090 
♦IUNI01 00 

* I UN 101 10 
* IUN I 01 20 

* 1 UN! 0 1 30 
* IUN I 01 4 0 
+ I UN I 0 1 50 
* IUN 101 60 
♦IUN 101 70 
♦ I UN I 01 80 

* I UN I 0 1 90 
* IUN 10200 
♦ I UN 10210 
♦ 1 UN I 022 0 
♦1 UN 10230 
♦ I UN I 0240 

* \ UN I 0250 
♦I UN I 0260 
♦ IUN I 0270 
♦ IUN 10280 
♦l UN 10290 

* I UN 10300 
♦IUN I 031 0 

* IUN I 0320 
* IUN I 0330 
♦1 UN I 0340 

* IUN I 0350 
♦ IUN 10360 

* I UN 1 0370 

* I UN 1 C380 
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c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

.c* 

c* 

C* ■ 

c* 

c* 

.c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* . 

c* 

c* 

C* 

c* 

c* 

c* 

c* 

c* 

c* 

c* ' 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c**-** 


c 


I ORDER 


xo 


1PT 


IE«R 


INTERPOLATION PARAMETER SUPPLIED BY THE USER. * 

* 

= 0 ZERO ORDER I NTERPOL AT I ON9 THE FIRST FUNCTION * 

VALUE IN EACH DEPENDENT VARIABLE TABLE IS * 

ASSIGNED TO THE CORRESPONDING MEMBER OF THE YO * 

ARRAY. the functional value is estimated to * 

REMAIN CONSTANT AND EQUAL TO THE NEAREST KNOWN • * 

FUNCTION VALUE. * 

* 

THE INPUT POINT AT WHICH INTERPOLATION WILL BE * 

PERFORMED. * 

* 

A ONE-DIMENSIONAL ARRAY DIMENSIONED (NTAB) IN THE * 
CALLING PEOGRAM. UPON RETURN THE ARRAY CONTAINS THE * 
ESTIMATED VALUE OF EACH FUNCTION AT XO* * 

* 

ON THE FIRST CALL 1PT MUST BE INITIALIZED TO -1 SO * 
THAT MONOTONICITY WILL BE CHECKED. UPON LEAVING THE * 
ROUTINE IPT EQUALS THE VALUE OF THE INDEX OF THE X * 

value preceding XO unless extrapolation WAS * 

PERFORMED* IN THAT CASE THE VALUE OF IPT IS * 

RETURNED AS9 * 

sU DENOTES XO *LT. X ( ! ) IF THE X ARRAY IS IN * 

INCREASING ORDER AND X(l) .GT. XO IF THE X ARRAY * 
IS IN DECREASING ORDER* * 

=N DENOTES XO *GT • X(N) IF THE X ARRAY IS IN * 

INCREASING ORDER AND XO .LT. X(N> IF THE X ARRAY * 

IS IN DECREASING ORDER. * 

* 

ON SUBSEQUENT CALLS* IPT IS USED AS A POINTER TO * 

BEGIN THE SEARCH FOR XO. * 

* 

ERROR PARAMETER GENERATED BY THE ROUTINE * 

=0 NORMAL return * 

=J THE J-TH ELEMENT OF THE X ARRAY IS OUT OF ORDER * 
=-l ZERO ORDER INTERPOLATION PERFORMED BECAUSE * 

I ORDER =0. * 

=-2 ZERO ORDER INTERPOLATION PERFORMED BECAUSE ONLY * 

ONE POINT WAS IN X ARRAY. * 

s-3. NO INTERPOLAT ION U/AS PERFORMED BECAUSE * 

INSUFFICIENT POINTS WERE SUPPLIED FOR SECOND * 

ORDER INTERPOLATION. * 

=-4 EXTRAPOLATION WAS PERFORMED * 

* 

UPON RETURN the PARAMETER IERR SHOULD BE TESTED IN * 
The CALLING PROGRAM. * 


REQUIRED ROUTINES 
SOURCE 

LANGUAGE 

DATE RELEASED 
LATEST REVISION 


NONE * 

* 

CMPB ROUTINE MTLUP MODIFIED * 

BY COMPUTER SCIENCES CORPORATION* 

* 

FORTRAN * 

* 

* 

AUGUST 1 « 1 973 * 

* 

JUNE 9* 1975 * 


DIMENSION X ( 1 ) • Y (NMAX* 1 ) . YO ( I ) 

’ NM1 =N- 1 
ifrr=o 
J= 1 

TEST FOR ZERO ORDER i NTERPOL A T I ON 

DELX = X ( 2 )-X ( 1 ) 

IF < 1 ORDER .EQ. O) GO TO 1U 
IF (N.LT. 2) GO TO 20 
GO To 50 


UN 10390 
UN I 0400 
UN I 04 1 0 
UN I 0420 
UN I 0430 
UN I 0440 
UN I 0450 
UN I 0460 
UN 104 70 
UN I 0480 
UN 1 0490 
UN I 0500 
UN I 051 0 
UN I 0520 
UN I 0530 
UN I 0540 
UN 1 0550 
UN I 0560 
UN I 0570 
UN I 0580 
UN I 0590 
UN I 0600 
UN I 06 1 0 
UN 1 0620 
UN I 0630 
UN I 0640 
UN 1 06^0 
UN I 0660 
UN 10670 
UN I 0680 
UN 10690 
UN 10700 
UN 1071 0 
UN 1 0720 
UN I 0730 
UN I 0740 
UN I 0750 
UN I 0760 
UN 10770 
UN I 0780 
UN 10790 
UN IOBOO 
UN I 08 10 
UN I0B20 
UN I 0830 
UN I 0840 
UN I 0850 
UN I 0860 
UN I 0870 
UN I 0880 
UN I 0890 
UN 1 0900 
UN I 09 1 0 
UN I 0920 
UN I 0930 
UN I 094 0 
UN I 0950 
UN 1 0960 
UN I 0970 
UN I 0980 
UN I 0990 
UN I 1 000 
UNI into 
UN I 1 020 
UN I 1 030 
UN l 1 040 
UN I 1 0B0 
UN I 1 060 
UNI 1 070 
UN I 1-080 
UNI 1 090 
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10 IE«R=-1 1UNI 

CO TO TO IUMI 

20 !ERR=-2 IUNI 

30 DO 40 NT=1.NTAB IUNl 

YOfNT )=Y ( 1 .NT ) 1 UN I 

40 CONTINUE 1UNI 

Return iuni 

50 IF (IPT .GT. -1) GO TO 65 I UN I 

0 IUNl 

C CHECK FOR TABLE OF NODE POINTS BEING STRICTLY MONOTONIC IUNl 

C THE SIGN OF OELX SIGNIFIES WHETHER TABLE IS IN IUNl 

C INCREASING OR DECREASING ORDER. IUNl 

c IUNl 

IF (OELX .EO. 0) GO TO 190 IUNl 

IF (N .EO. 2) GO TO 65 IUNl 

c IUNl 

C check FOR SIGN CONSISTENCY IN THE DIFFERENCES OF IUNl 

SUBSEQUENT PAIRS IUNl 

IUNl 

DO 60 J = 2 . NM 1 [UNI 

IF (OELX * ( X ( JF l ) -X ( J ) ) ) 190.190.60 IUNl 

60 CONTINUE IUNl 

I UN I 

ipt is initialized to dE within the interval iuni 

I UN I 

65 IF (IPT .LT. 1) I PT = I IUNI 

IF (IPT .GT. MM 1 1 1PT=NM1 IUNl 

IN= SIGN ( I .0 .OELX *< XU-X(IPT)l) IUNI 

70 P= XIIPTI - *0 IUNI 

IF (P* ( X ( IPT +1)- XO ) ) 90.160. 80 IUNI 

80 IPT = I P T +IN IUNI 

1 UN I 

TEST TO SEE IF IT IS NECCESARY TO EXTRAPOLATE IUNI 

IUNI 

IF (IPT.GT.O .ANO. IPT .(_T. N) GO TO 70 IUNl 

IERR=-4 IUNI 

I PT= IPT- IN IUNI 

c IUNI 

C TEST FOR ORDER OF INTERPOLATION IUNI 

C IUNI 

C IUNl 

90 IF (I ORDER . GT. 1 ) GO TO 120 IUNI 

C IUNI 

c first order INTERPOLATION IUNI 

c IUNI 

I PT I = I PT+ 1 I UN I 

XTMP 1 =X0-X ( IPT) IUNI 

XTMP2=X ( IPT1 )-X ( IPT ) IUNl 

XTMP1 =XTMPl /XTMP2 1 UN I 

DO 100 NT=l.NTAe) IUNI 

YTMP=Y< IPT1 .NT )-Y( IPT, NT ) IUNl 

YC(NT)=YI IPT.NT |+YTMP*XTMP1 IUNl 

1 00 CONTINUE IUNI 

IF ( I ERR .EO. -4) IPT=|PT+IN IUNI 

RETURN IUNI 

C I UN I 

C SECOND ORDER INTERPOLATION IUNl 

C IUNl 

120 IF (N .EO. 2) GO TO 200 IUNl 

C IUNl 

C CHOOSING A THIRD POINT SO AS TO MINIMIZE THE DISTANCE IUNl 

C BETWEEN THE THREE POINTS USED TO INTERPOLATE IUNl 

C IUNI 

IF (IPT .EO. NM1 ) GO TO 140 IUNl 

IF (IPT .EO. 1) GO TO 130 IUNl 

A1 =ABS(X0-X( IPT— 1 ) ) IUNl 

A2=A3S(X( IPT+21-X0) IUNI 

IF ( A 1 -A2 ) 140. 1 2U. 1 30 IUNI 

130 L= IPT IUNl 

GO TO 150 IUNl 

1 4 O L = I P T - 1 IUNl 

ISO V 1 -X (L ) — XO IUNl 


1 I 00 
1110 
1 120 
1 1 30 
1 140 
1 150 
1 I 60 
1 1 70 
1 180 
1 190 
1200 
1210 
1220 
1230 
1240 
1250 
1 260 
1 270 
1280 
1 290 
1 300 
1310 
1 320 
1 330 
1340 
1 350 
1 360 
1 370 
1380 
1 390 
1400 
14 10 
1420 
1430 
1440 
1 450 
1460 
1470 
1480 
1490 
1500 
1510 
1520 
1 530 
1540 
I 550 
1 560 
1 570 
1580 
1590 
1600 
1610 
1620 
1 630 
1640 
1 650 
1660 
1670 
1680 
1690 
1 700 
1710 
1 720 
I 730 
1 740 
1 750 
1 760 
1 770 
1 780 
1 790 
1800 
1810 
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V2=X (L+l )-X0 
V3=X(L+2)-X0 
00 I 60 NT = I •NT/SB 

YY1=(Y(L.NT) * V2 - Y(L+1.NT) * Vl)/(X(L+1) - X(L)) 

YY2= <Y(L+1 iNT ) #V3— Y <L+2*NT ) *V2 ) / ( X < L+2 1 - X (L + 1 )) 

Y0 (NT >= <YY1 *V3-YY2*V1 ) / ( X < L +2 ) -X ( L ) ) 

260 CONTINUE 

IF (I ERR *EQ. -4) I PT= I PT + IN 
RETURN 

1 80 . I F ( P .NE* O) I PT = I P T +1 

DO 1 35 NT = 1 • NT AS 

YCMNT )=Y< IPT*NT ) 

185 CONTINUE 

RETURN 
C 

IERR IS SET TO THE SUBSCRIPT OF THE MEMBFR OF Thf TABLE 
WHICH IS OUT OF ORDER 
C 

190 IERR=J +1 
RETURN 

200 I ERR=-3 
RETURN 
END 


I UN I 1820 
I UN I 1830 
I UN I 1840 
I UN I 1850 
I UN I I860 
I UN I 1870 
I UNI ! 880 
I UN I 1890 
I UN I 1900 
I UN I 1 91 0 
! UN I 1 920 
I UN I 1 930 
I UN I 1940 
I UN I 1950 
I UN I 1 960 
I UN I 1 970 
1 UN 1 1960 
I UN I 1 990 
I UNI 2000 
I UN I 20 1 0 
I UN I 2020 
I UN 12030 
I UN I 2040 
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TABLE 2.- BUILDUP-FACTOR PARAMETERS 


E, GeV 

Buildup -factor parameters 
for dose equivalent 

Buildup -factor parameters 
for absorbed dose 

a l 

a 2 

a 3 

a 4 

a l 

a 2 

a 3 

a 4 

0.03 

*1.00 

*0 

*0 

*0 

*1.00 

*0 

*0 


.06 

*1.20 

*0 

*0 

*.0130 

*1.07 

*.010 

*0 

*.010 

.10 

1.40 

.020 

0 

.0300 

1.10 

.040 

0 

.026 

.15 

*1.50 

*.070 

*0 

*.0385 

*1.12 

*.060 

*0 

*.031 

.20 

1.60 

.090 

0 

.0400 

1.15 

.062 

0 

.032 

.30 

1.70 

.110 

0 

.0330 

1.20 

.068 

0 

.026 

.40 

1.90 

.130 

0 

.0228 

1.24 

.071 

0 

.0228 

.73 

3.40 

.156 

.00035 

.0150 

1.40 

.090 


.0150 

1.2 

*4.32 

*.167 

*.00145 

*.0130 

*1.67 

*.094 


*.0122 

1.5 

4.60 

.170 

.00250 

.0120 

1.80 

.095 


.0120 

3.0 

5.35 

.190 

.00300 

.0100 

2.00 

.100 



10.0 

6.20 

.280 

.00350 

.0100 

2.30 

.111 

.00205 



^Values obtained by interpolation. 
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Figure 3.- Contribution of nuclear reactions to absorbed Figure 4.- Contribution of nuclear reactions t 
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